function [T0UU T1UU T2UU] = h2stressenergy(params, u, udt, udtdt)

    Nx = params.Nx;
    Ny = params.Ny; 
    Nt = Nx.*Ny; 

    T  = reshape(u(0*Nt+1 : 1*Nt), Nx, Ny);
    ux = reshape(u(1*Nt+1 : 2*Nt), Nx, Ny);
    uy = reshape(u(2*Nt+1 : 3*Nt), Nx, Ny);
    ut = reshape(u(3*Nt+1 : 4*Nt), Nx, Ny);

    uU = zeros(Nx, Ny, 3);
    uU(:, :, 1) = ut;
    uU(:, :, 2) = ux;
    uU(:, :, 3) = uy;


    Tdt  = reshape(udt(0*Nt+1 : 1*Nt), Nx, Ny);
    uxdt = reshape(udt(1*Nt+1 : 2*Nt), Nx, Ny);
    uydt = reshape(udt(2*Nt+1 : 3*Nt), Nx, Ny);
    utdt = reshape(udt(3*Nt+1 : 4*Nt), Nx, Ny);

    uxdtdt = reshape(udtdt(1*Nt+1 : 2*Nt), Nx, Ny);
    uydtdt = reshape(udtdt(2*Nt+1 : 3*Nt), Nx, Ny);
    utdtdt = reshape(udtdt(3*Nt+1 : 4*Nt), Nx, Ny);

    Tdx = dx(T);
    Tdy = dy(T);

    utdx = dx(ut);
    utdy = dy(ut);

    uxdx = dx(ux);
    uxdy = dy(ux);

    uydx = dx(uy);
    uydy = dy(uy);


    utdtdx  = dx(utdt);
    utdtdy  = dy(utdt);
    utdxdx  = dx(utdx);
    utdxdy  = dy(utdx);
    utdydy  = dy(utdy);

    uxdtdx  = dx(uxdt);
    uxdtdy  = dy(uxdt);
    uxdxdx  = dx(uxdx);
    uxdxdy  = dy(uxdx);
    uxdydy  = dy(uxdy);

    uydtdx  = dx(uydt);
    uydtdy  = dy(uydt);
    uydxdx  = dx(uydx);
    uydxdy  = dy(uydx);
    uydydy  = dy(uydy);

    sUU = zeros(params.Nx, params.Ny, 3, 3);  % sigmaUU
    sUU(:,:,1,1) = (1/2).*((-1).*((-1)+ut.^2).^2.*utdt+uxdx+ux.^2.*uxdx+ux.*uxdy.*uy+ut.*(utdx.*ux+(-1).*ux.*uxdt+uy.*(utdy+(-1).*uydt))+ut.^3.*((-1).*utdx.*ux+ux.*uxdt+uy.*((-1).*utdy+uydt))+ux.*uy.*uydx+uydy+uy.^2.*uydy+ut.^2.*(((-1)+ux.^2).*uxdx+ux.*uy.*(uxdy+uydx)+((-1)+uy.^2).*uydy));
    sUU(:,:,1,2) = (1/2).*((-1).*ut.^3.*utdt.*ux+(-1).*((-1)+ut.^2).*utdx.*(1+ux.^2)+(-1).*uxdt+(-1).*ux.^2.*uxdt+utdy.*ux.*uy+(-1).*ux.*uy.*uydt+ut.^2.*((1+ux.^2).*uxdt+ux.*uy.*((-1).*utdy+uydt))+ut.*(utdt.*ux+ux.^3.*uxdx+uy.*(uxdy+uydx)+ux.^2.*uy.*(uxdy+uydx)+ux.*(uxdx+((-1)+uy.^2).*uydy)));
    sUU(:,:,1,3) = (1/2).*((-1).*ut.^3.*utdt.*uy+utdx.*ux.*uy+(-1).*ux.*uxdt.*uy+(-1).*((-1)+ut.^2).*utdy.*(1+uy.^2)+(-1).*uydt+(-1).*uy.^2.*uydt+ut.^2.*((-1).*utdx.*ux.*uy+ux.*uxdt.*uy+uydt+uy.^2.*uydt)+ut.*(ux.^2.*uxdx.*uy+ux.*(1+uy.^2).*(uxdy+uydx)+uy.*(utdt+(-1).*uxdx+uydy+uy.^2.*uydy)));

    sUU(:,:,2,2) = (1/2).*((-1).*utdt.*(1+ux.^2+ut.^2.*((-1)+ux.^2))+uxdx+2.*ux.^2.*uxdx+ux.^4.*uxdx+ux.*uxdy.*uy+ux.^3.*uxdy.*uy+ut.*((-1).*utdx.*(ux+ux.^3)+ux.*uxdt+ux.^3.*uxdt+uy.*(utdy+(-1).*uydt)+ux.^2.*uy.*((-1).*utdy+uydt))+ux.*uy.*uydx+ux.^3.*uy.*uydx+(-1).*uydy+(-1).*ux.^2.*uydy+(-1).*uy.^2.*uydy+ux.^2.*uy.^2.*uydy);
    sUU(:,:,2,3) = (1/2).*((-1).*utdt.*ux.*uy+(-1).*ut.^2.*utdt.*ux.*uy+ux.*uxdx.*uy+ux.^3.*uxdx.*uy+(1+ux.^2).*uxdy.*(1+uy.^2)+ut.*((-1).*utdx.*(1+ux.^2).*uy+uxdt.*uy+ux.^2.*uxdt.*uy+(-1).*utdy.*ux.*(1+uy.^2)+ux.*uydt+ux.*uy.^2.*uydt)+uydx+ux.^2.*uydx+uy.^2.*uydx+ux.^2.*uy.^2.*uydx+ux.*uy.*uydy+ux.*uy.^3.*uydy);

    sUU(:,:,3,3) = (1/2).*((-1).*uxdx+(-1).*ux.^2.*uxdx+ux.*uxdy.*uy+(-1).*uxdx.*uy.^2+ux.^2.*uxdx.*uy.^2+ux.*uxdy.*uy.^3+(-1).*utdt.*(1+uy.^2+ut.^2.*((-1)+uy.^2))+ut.*(ux.*uxdt.*((-1)+uy.^2)+utdx.*(ux+(-1).*ux.*uy.^2)+(-1).*uy.*(1+uy.^2).*(utdy+(-1).*uydt))+ux.*uy.*uydx+ux.*uy.^3.*uydx+uydy+2.*uy.^2.*uydy+uy.^4.*uydy);


    sUUd = zeros(params.Nx, params.Ny, 3, 3, 3);  % sigmaUU,d
    sUUd(:,:,1,1,1) = (1/2).*((-1).*utdtdt+(-1).*ut.^4.*utdtdt+utdt.*utdx.*ux+(-1).*utdt.*ux.*uxdt+uxdtdx+ux.^2.*uxdtdx+2.*ux.*uxdt.*uxdx+utdt.*utdy.*uy+ux.*uxdtdy.*uy+uxdt.*uxdy.*uy+ux.*uxdy.*uydt+(-1).*utdt.*uy.*uydt+ut.^3.*((-4).*utdt.^2+(-1).*utdtdx.*ux+(-1).*utdx.*uxdt+uxdt.^2+ux.*uxdtdt+(-1).*utdtdy.*uy+(-1).*utdy.*uydt+uydt.^2+uy.*uydtdt)+ux.*uy.*uydtdx+uydtdy+uy.^2.*uydtdy+uxdt.*uy.*uydx+ux.*uydt.*uydx+2.*uy.*uydt.*uydy+ut.^2.*(2.*utdtdt+(-1).*uxdtdx+ux.^2.*uxdtdx+2.*ux.*uxdt.*uxdx+ux.*uxdtdy.*uy+uxdt.*uxdy.*uy+ux.*uxdy.*uydt+(-3).*utdt.*(utdx.*ux+(-1).*ux.*uxdt+utdy.*uy+(-1).*uy.*uydt)+ux.*uy.*uydtdx+(-1).*uydtdy+uy.^2.*uydtdy+uxdt.*uy.*uydx+ux.*uydt.*uydx+2.*uy.*uydt.*uydy)+ut.*(4.*utdt.^2+utdtdx.*ux+utdx.*uxdt+(-1).*uxdt.^2+(-1).*ux.*uxdtdt+utdtdy.*uy+utdy.*uydt+(-1).*uydt.^2+(-1).*uy.*uydtdt+2.*utdt.*(((-1)+ux.^2).*uxdx+ux.*uy.*(uxdy+uydx)+((-1)+uy.^2).*uydy)));
    sUUd(:,:,1,1,2) = (1/2).*((-1).*((-1)+ut.^2).^2.*utdtdx+utdx.^2.*ux+(-1).*utdx.*ux.*uxdt+2.*ux.*uxdx.^2+uxdxdx+ux.^2.*uxdxdx+utdx.*utdy.*uy+ux.*uxdxdy.*uy+uxdx.*uxdy.*uy+(-1).*utdx.*uy.*uydt+ux.*uxdy.*uydx+uxdx.*uy.*uydx+ux.*uydx.^2+ut.^3.*((-4).*utdt.*utdx+(-1).*utdxdx.*ux+ux.*uxdtdx+(-1).*utdx.*uxdx+uxdt.*uxdx+(-1).*utdxdy.*uy+uy.*uydtdx+(-1).*utdy.*uydx+uydt.*uydx)+ux.*uy.*uydxdx+uydxdy+uy.^2.*uydxdy+2.*uy.*uydx.*uydy+ut.*(4.*utdt.*utdx+utdxdx.*ux+(-1).*ux.*uxdtdx+(-1).*utdx.*uxdx+2.*utdx.*ux.^2.*uxdx+(-1).*uxdt.*uxdx+utdxdy.*uy+2.*utdx.*ux.*uxdy.*uy+(-1).*uy.*uydtdx+utdy.*uydx+2.*utdx.*ux.*uy.*uydx+(-1).*uydt.*uydx+(-2).*utdx.*uydy+2.*utdx.*uy.^2.*uydy)+ut.^2.*((-3).*utdx.^2.*ux+(-1).*uxdxdx+ux.^2.*uxdxdx+uxdx.*uxdy.*uy+3.*utdx.*(ux.*uxdt+uy.*((-1).*utdy+uydt))+uxdx.*uy.*uydx+ux.*(2.*uxdx.^2+uxdxdy.*uy+uxdy.*uydx+uydx.^2+uy.*uydxdx)+(-1).*uydxdy+uy.^2.*uydxdy+2.*uy.*uydx.*uydy));
    sUUd(:,:,1,1,3) = (1/2).*((-1).*((-1)+ut.^2).^2.*utdtdy+utdx.*utdy.*ux+(-1).*utdy.*ux.*uxdt+uxdxdy+ux.^2.*uxdxdy+2.*ux.*uxdx.*uxdy+utdy.^2.*uy+uxdy.^2.*uy+ux.*uxdydy.*uy+(-1).*utdy.*uy.*uydt+uxdy.*uy.*uydx+ux.*uy.*uydxdy+ux.*uxdy.*uydy+ux.*uydx.*uydy+2.*uy.*uydy.^2+ut.*(4.*utdt.*utdy+utdxdy.*ux+(-1).*ux.*uxdtdy+(-2).*utdy.*uxdx+2.*utdy.*ux.^2.*uxdx+utdx.*uxdy+(-1).*uxdt.*uxdy+utdydy.*uy+2.*utdy.*ux.*uxdy.*uy+(-1).*uy.*uydtdy+2.*utdy.*ux.*uy.*uydx+(-1).*utdy.*uydy+2.*utdy.*uy.^2.*uydy+(-1).*uydt.*uydy)+ut.^3.*((-4).*utdt.*utdy+(-1).*utdxdy.*ux+ux.*uxdtdy+(-1).*utdx.*uxdy+uxdt.*uxdy+(-1).*utdydy.*uy+uy.*uydtdy+(-1).*utdy.*uydy+uydt.*uydy)+uydydy+uy.^2.*uydydy+ut.^2.*((-3).*utdx.*utdy.*ux+3.*utdy.*ux.*uxdt+(-1).*uxdxdy+ux.^2.*uxdxdy+2.*ux.*uxdx.*uxdy+(-3).*utdy.^2.*uy+uxdy.^2.*uy+ux.*uxdydy.*uy+3.*utdy.*uy.*uydt+uxdy.*uy.*uydx+ux.*uy.*uydxdy+ux.*uxdy.*uydy+ux.*uydx.*uydy+2.*uy.*uydy.^2+(-1).*uydydy+uy.^2.*uydydy));

    sUUd(:,:,1,2,1) = (1/2).*(utdt.^2.*ux+(-1).*((-1)+ut.^2).*utdtdx.*(1+ux.^2)+2.*utdx.*ux.*uxdt+(-2).*ux.*uxdt.^2+(-1).*ut.^3.*(utdtdt.*ux+utdt.*uxdt)+(-1).*uxdtdt+(-1).*ux.^2.*uxdtdt+utdt.*ux.*uxdx+utdt.*ux.^3.*uxdx+utdtdy.*ux.*uy+utdy.*uxdt.*uy+utdt.*uxdy.*uy+utdt.*ux.^2.*uxdy.*uy+utdy.*ux.*uydt+(-1).*uxdt.*uy.*uydt+(-1).*ux.*uydt.^2+(-1).*ux.*uy.*uydtdt+ut.^2.*((-3).*utdt.^2.*ux+(-2).*utdx.*ux.*uxdt+2.*ux.*uxdt.^2+uxdtdt+ux.^2.*uxdtdt+(-1).*utdtdy.*ux.*uy+(-1).*utdy.*uxdt.*uy+(-1).*utdy.*ux.*uydt+uxdt.*uy.*uydt+ux.*uydt.^2+ux.*uy.*uydtdt)+utdt.*uy.*uydx+utdt.*ux.^2.*uy.*uydx+(-1).*utdt.*ux.*uydy+utdt.*ux.*uy.^2.*uydy+ut.*(utdtdt.*ux+ux.*uxdtdx+ux.^3.*uxdtdx+uxdt.*uxdx+3.*ux.^2.*uxdt.*uxdx+uxdtdy.*uy+ux.^2.*uxdtdy.*uy+2.*ux.*uxdt.*uxdy.*uy+uxdy.*uydt+ux.^2.*uxdy.*uydt+utdt.*((-2).*utdx.*(1+ux.^2)+(3+2.*ux.^2).*uxdt+2.*ux.*uy.*((-1).*utdy+uydt))+uy.*uydtdx+ux.^2.*uy.*uydtdx+(-1).*ux.*uydtdy+ux.*uy.^2.*uydtdy+2.*ux.*uxdt.*uy.*uydx+uydt.*uydx+ux.^2.*uydt.*uydx+(-1).*uxdt.*uydy+uxdt.*uy.^2.*uydy+2.*ux.*uy.*uydt.*uydy));
    sUUd(:,:,1,2,2) = (1/2).*(utdxdx+utdt.*utdx.*ux+utdxdx.*ux.^2+(-1).*uxdtdx+(-1).*ux.^2.*uxdtdx+3.*utdx.*ux.*uxdx+utdx.*ux.^3.*uxdx+(-2).*ux.*uxdt.*uxdx+(-1).*ut.^3.*(utdtdx.*ux+utdt.*uxdx)+utdxdy.*ux.*uy+utdy.*uxdx.*uy+utdx.*uxdy.*uy+utdx.*ux.^2.*uxdy.*uy+(-1).*uxdx.*uy.*uydt+(-1).*ux.*uy.*uydtdx+utdy.*ux.*uydx+utdx.*uy.*uydx+utdx.*ux.^2.*uy.*uydx+(-1).*ux.*uydt.*uydx+ut.^2.*((-3).*utdt.*utdx.*ux+(-1).*utdxdx.*(1+ux.^2)+uxdtdx+ux.^2.*uxdtdx+(-2).*utdx.*ux.*uxdx+2.*ux.*uxdt.*uxdx+(-1).*utdxdy.*ux.*uy+(-1).*utdy.*uxdx.*uy+uxdx.*uy.*uydt+ux.*uy.*uydtdx+(-1).*utdy.*ux.*uydx+ux.*uydt.*uydx)+(-1).*utdx.*ux.*uydy+utdx.*ux.*uy.^2.*uydy+ut.*(utdtdx.*ux+(-2).*utdx.^2.*(1+ux.^2)+utdt.*uxdx+uxdx.^2+3.*ux.^2.*uxdx.^2+ux.*uxdxdx+ux.^3.*uxdxdx+uxdxdy.*uy+ux.^2.*uxdxdy.*uy+2.*ux.*uxdx.*uxdy.*uy+2.*utdx.*(uxdt+ux.^2.*uxdt+(-1).*utdy.*ux.*uy+ux.*uy.*uydt)+uxdy.*uydx+ux.^2.*uxdy.*uydx+2.*ux.*uxdx.*uy.*uydx+uydx.^2+ux.^2.*uydx.^2+uy.*uydxdx+ux.^2.*uy.*uydxdx+(-1).*ux.*uydxdy+ux.*uy.^2.*uydxdy+(-1).*uxdx.*uydy+uxdx.*uy.^2.*uydy+2.*ux.*uy.*uydx.*uydy));
    sUUd(:,:,1,2,3) = (1/2).*(utdt.*utdy.*ux+(-1).*((-1)+ut.^2).*utdxdy.*(1+ux.^2)+(-1).*uxdtdy+(-1).*ux.^2.*uxdtdy+utdy.*ux.*uxdx+utdy.*ux.^3.*uxdx+2.*utdx.*ux.*uxdy+(-2).*ux.*uxdt.*uxdy+(-1).*ut.^3.*(utdtdy.*ux+utdt.*uxdy)+utdydy.*ux.*uy+2.*utdy.*uxdy.*uy+utdy.*ux.^2.*uxdy.*uy+(-1).*uxdy.*uy.*uydt+(-1).*ux.*uy.*uydtdy+utdy.*uy.*uydx+utdy.*ux.^2.*uy.*uydx+utdy.*ux.*uy.^2.*uydy+(-1).*ux.*uydt.*uydy+ut.^2.*((-3).*utdt.*utdy.*ux+uxdtdy+ux.^2.*uxdtdy+(-2).*utdx.*ux.*uxdy+2.*ux.*uxdt.*uxdy+(-1).*utdydy.*ux.*uy+(-1).*utdy.*uxdy.*uy+uxdy.*uy.*uydt+ux.*uy.*uydtdy+(-1).*utdy.*ux.*uydy+ux.*uydt.*uydy)+ut.*(utdtdy.*ux+(-2).*utdx.*utdy.*(1+ux.^2)+2.*utdy.*uxdt+2.*utdy.*ux.^2.*uxdt+ux.*uxdxdy+ux.^3.*uxdxdy+utdt.*uxdy+uxdx.*uxdy+3.*ux.^2.*uxdx.*uxdy+(-2).*utdy.^2.*ux.*uy+2.*ux.*uxdy.^2.*uy+uxdydy.*uy+ux.^2.*uxdydy.*uy+2.*utdy.*ux.*uy.*uydt+2.*ux.*uxdy.*uy.*uydx+uy.*uydxdy+ux.^2.*uy.*uydxdy+ux.^2.*uxdy.*uydy+uxdy.*uy.^2.*uydy+uydx.*uydy+ux.^2.*uydx.*uydy+2.*ux.*uy.*uydy.^2+(-1).*ux.*uydydy+ux.*uy.^2.*uydydy));

    sUUd(:,:,1,3,1) = (1/2).*(utdt.*ux.*uxdy+utdt.^2.*uy+utdtdx.*ux.*uy+utdx.*uxdt.*uy+(-1).*uxdt.^2.*uy+(-1).*ux.*uxdtdt.*uy+(-1).*utdt.*uxdx.*uy+utdt.*ux.^2.*uxdx.*uy+utdt.*ux.*uxdy.*uy.^2+(-1).*((-1)+ut.^2).*utdtdy.*(1+uy.^2)+utdx.*ux.*uydt+(-1).*ux.*uxdt.*uydt+2.*utdy.*uy.*uydt+(-2).*uy.*uydt.^2+(-1).*ut.^3.*(utdtdt.*uy+utdt.*uydt)+(-1).*uydtdt+(-1).*uy.^2.*uydtdt+ut.^2.*((-3).*utdt.^2.*uy+(-1).*utdtdx.*ux.*uy+(-1).*utdx.*uxdt.*uy+uxdt.^2.*uy+ux.*uxdtdt.*uy+(-1).*utdx.*ux.*uydt+ux.*uxdt.*uydt+(-2).*utdy.*uy.*uydt+2.*uy.*uydt.^2+uydtdt+uy.^2.*uydtdt)+utdt.*ux.*uydx+utdt.*ux.*uy.^2.*uydx+utdt.*uy.*uydy+utdt.*uy.^3.*uydy+ut.*(uxdt.*uxdy+utdtdt.*uy+(-1).*uxdtdx.*uy+uxdt.*uxdy.*uy.^2+(-1).*uxdx.*uydt+ux.^2.*(uxdtdx.*uy+uxdx.*uydt)+utdt.*((-2).*utdx.*ux.*uy+2.*ux.*uxdt.*uy+(-2).*utdy.*(1+uy.^2)+3.*uydt+2.*uy.^2.*uydt)+uy.*uydtdy+uy.^3.*uydtdy+uxdt.*uydx+uxdt.*uy.^2.*uydx+ux.*(uxdtdy+2.*uxdt.*uxdx.*uy+uxdtdy.*uy.^2+2.*uxdy.*uy.*uydt+uydtdx+uy.^2.*uydtdx+2.*uy.*uydt.*uydx)+uydt.*uydy+3.*uy.^2.*uydt.*uydy));
    sUUd(:,:,1,3,2) = (1/2).*(utdx.*ux.*uxdy+utdt.*utdx.*uy+utdxdx.*ux.*uy+(-1).*ux.*uxdtdx.*uy+utdx.*ux.^2.*uxdx.*uy+(-1).*uxdt.*uxdx.*uy+utdx.*ux.*uxdy.*uy.^2+(-1).*((-1)+ut.^2).*utdxdy.*(1+uy.^2)+(-1).*uydtdx+(-1).*uy.^2.*uydtdx+2.*utdx.*ux.*uydx+(-1).*ux.*uxdt.*uydx+2.*utdy.*uy.*uydx+utdx.*ux.*uy.^2.*uydx+(-2).*uy.*uydt.*uydx+(-1).*ut.^3.*(utdtdx.*uy+utdt.*uydx)+ut.^2.*((-3).*utdt.*utdx.*uy+(-1).*utdxdx.*ux.*uy+ux.*uxdtdx.*uy+(-1).*utdx.*uxdx.*uy+uxdt.*uxdx.*uy+uydtdx+uy.^2.*uydtdx+(-1).*utdx.*ux.*uydx+ux.*uxdt.*uydx+(-2).*utdy.*uy.*uydx+2.*uy.*uydt.*uydx)+utdx.*uy.*uydy+utdx.*uy.^3.*uydy+ut.*(uxdx.*uxdy+utdtdx.*uy+(-2).*utdx.^2.*ux.*uy+(-1).*uxdxdx.*uy+uxdx.*uxdy.*uy.^2+2.*utdx.*(ux.*uxdt.*uy+(-1).*utdy.*(1+uy.^2)+uydt+uy.^2.*uydt)+utdt.*uydx+uxdx.*uy.^2.*uydx+ux.^2.*(uxdxdx.*uy+uxdx.*uydx)+ux.*(uxdxdy+2.*uxdx.^2.*uy+uxdxdy.*uy.^2+2.*uxdy.*uy.*uydx+2.*uy.*uydx.^2+uydxdx+uy.^2.*uydxdx)+uy.*uydxdy+uy.^3.*uydxdy+uydx.*uydy+3.*uy.^2.*uydx.*uydy));
    sUUd(:,:,1,3,3) = (1/2).*(utdydy+utdy.*ux.*uxdy+utdt.*utdy.*uy+utdxdy.*ux.*uy+(-1).*ux.*uxdtdy.*uy+(-1).*utdy.*uxdx.*uy+utdy.*ux.^2.*uxdx.*uy+utdx.*uxdy.*uy+(-1).*uxdt.*uxdy.*uy+utdydy.*uy.^2+utdy.*ux.*uxdy.*uy.^2+(-1).*uydtdy+(-1).*uy.^2.*uydtdy+utdy.*ux.*uydx+utdy.*ux.*uy.^2.*uydx+utdx.*ux.*uydy+(-1).*ux.*uxdt.*uydy+3.*utdy.*uy.*uydy+utdy.*uy.^3.*uydy+(-2).*uy.*uydt.*uydy+(-1).*ut.^3.*(utdtdy.*uy+utdt.*uydy)+ut.^2.*((-3).*utdt.*utdy.*uy+(-1).*utdxdy.*ux.*uy+ux.*uxdtdy.*uy+(-1).*utdx.*uxdy.*uy+uxdt.*uxdy.*uy+(-1).*utdydy.*(1+uy.^2)+uydtdy+uy.^2.*uydtdy+(-1).*utdx.*ux.*uydy+ux.*uxdt.*uydy+(-2).*utdy.*uy.*uydy+2.*uy.*uydt.*uydy)+ut.*(ux.*uxdydy+utdtdy.*uy+(-1).*uxdxdy.*uy+ux.^2.*uxdxdy.*uy+ux.*uxdydy.*uy.^2+(-2).*utdy.^2.*(1+uy.^2)+uxdy.^2.*(1+uy.^2)+2.*utdy.*((-1).*utdx.*ux.*uy+ux.*uxdt.*uy+uydt+uy.^2.*uydt)+ux.*uydxdy+ux.*uy.^2.*uydxdy+utdt.*uydy+(-1).*uxdx.*uydy+ux.^2.*uxdx.*uydy+2.*ux.*uy.*uydx.*uydy+uydy.^2+3.*uy.^2.*uydy.^2+uxdy.*((1+uy.^2).*uydx+2.*ux.*uy.*(uxdx+uydy))+uy.*uydydy+uy.^3.*uydydy));

    sUUd(:,:,2,2,1) = (1/2).*((-1).*utdt.*utdx.*ux+(-1).*utdt.*utdx.*ux.^3+(-1).*utdtdt.*(1+ux.^2)+(-1).*utdt.*ux.*uxdt+utdt.*ux.^3.*uxdt+ut.^2.*(utdtdt+(-1).*utdtdt.*ux.^2+(-2).*utdt.*ux.*uxdt)+uxdtdx+2.*ux.^2.*uxdtdx+ux.^4.*uxdtdx+4.*ux.*uxdt.*uxdx+4.*ux.^3.*uxdt.*uxdx+utdt.*utdy.*uy+(-1).*utdt.*utdy.*ux.^2.*uy+ux.*uxdtdy.*uy+ux.^3.*uxdtdy.*uy+uxdt.*uxdy.*uy+3.*ux.^2.*uxdt.*uxdy.*uy+ux.*uxdy.*uydt+ux.^3.*uxdy.*uydt+(-1).*utdt.*uy.*uydt+utdt.*ux.^2.*uy.*uydt+ut.*((-2).*utdt.^2.*((-1)+ux.^2)+(-1).*utdtdx.*(ux+ux.^3)+(-1).*utdx.*uxdt+(-3).*utdx.*ux.^2.*uxdt+uxdt.^2+3.*ux.^2.*uxdt.^2+ux.*uxdtdt+ux.^3.*uxdtdt+utdtdy.*uy+(-1).*utdtdy.*ux.^2.*uy+(-2).*utdy.*ux.*uxdt.*uy+utdy.*uydt+(-1).*utdy.*ux.^2.*uydt+2.*ux.*uxdt.*uy.*uydt+(-1).*uydt.^2+ux.^2.*uydt.^2+(-1).*uy.*uydtdt+ux.^2.*uy.*uydtdt)+ux.*uy.*uydtdx+ux.^3.*uy.*uydtdx+(-1).*uydtdy+(-1).*ux.^2.*uydtdy+(-1).*uy.^2.*uydtdy+ux.^2.*uy.^2.*uydtdy+uxdt.*uy.*uydx+3.*ux.^2.*uxdt.*uy.*uydx+ux.*uydt.*uydx+ux.^3.*uydt.*uydx+(-2).*ux.*uxdt.*uydy+2.*ux.*uxdt.*uy.^2.*uydy+(-2).*uy.*uydt.*uydy+2.*ux.^2.*uy.*uydt.*uydy);
    sUUd(:,:,2,2,2) = (1/2).*((-1).*utdx.^2.*ux+(-1).*utdx.^2.*ux.^3+(-1).*utdtdx.*(1+ux.^2+ut.^2.*((-1)+ux.^2))+utdx.*ux.*uxdt+utdx.*ux.^3.*uxdt+(-2).*utdt.*ux.*uxdx+(-2).*ut.^2.*utdt.*ux.*uxdx+4.*ux.*uxdx.^2+4.*ux.^3.*uxdx.^2+uxdxdx+2.*ux.^2.*uxdxdx+ux.^4.*uxdxdx+utdx.*utdy.*uy+(-1).*utdx.*utdy.*ux.^2.*uy+ux.*uxdxdy.*uy+ux.^3.*uxdxdy.*uy+uxdx.*uxdy.*uy+3.*ux.^2.*uxdx.*uxdy.*uy+(-1).*utdx.*uy.*uydt+utdx.*ux.^2.*uy.*uydt+ux.*uxdy.*uydx+ux.^3.*uxdy.*uydx+uxdx.*uy.*uydx+3.*ux.^2.*uxdx.*uy.*uydx+ux.*uydx.^2+ux.^3.*uydx.^2+ut.*((-2).*utdt.*utdx.*((-1)+ux.^2)+(-1).*utdxdx.*(ux+ux.^3)+ux.*uxdtdx+ux.^3.*uxdtdx+(-1).*utdx.*uxdx+(-3).*utdx.*ux.^2.*uxdx+uxdt.*uxdx+3.*ux.^2.*uxdt.*uxdx+utdxdy.*uy+(-1).*utdxdy.*ux.^2.*uy+(-2).*utdy.*ux.*uxdx.*uy+2.*ux.*uxdx.*uy.*uydt+(-1).*uy.*uydtdx+ux.^2.*uy.*uydtdx+utdy.*uydx+(-1).*utdy.*ux.^2.*uydx+(-1).*uydt.*uydx+ux.^2.*uydt.*uydx)+ux.*uy.*uydxdx+ux.^3.*uy.*uydxdx+(-1).*uydxdy+(-1).*ux.^2.*uydxdy+(-1).*uy.^2.*uydxdy+ux.^2.*uy.^2.*uydxdy+(-2).*ux.*uxdx.*uydy+2.*ux.*uxdx.*uy.^2.*uydy+(-2).*uy.*uydx.*uydy+2.*ux.^2.*uy.*uydx.*uydy);
    sUUd(:,:,2,2,3) = (1/2).*((-1).*utdx.*utdy.*ux+(-1).*utdx.*utdy.*ux.^3+(-1).*utdtdy.*(1+ux.^2+ut.^2.*((-1)+ux.^2))+utdy.*ux.*uxdt+utdy.*ux.^3.*uxdt+uxdxdy+2.*ux.^2.*uxdxdy+ux.^4.*uxdxdy+(-2).*utdt.*ux.*uxdy+(-2).*ut.^2.*utdt.*ux.*uxdy+4.*ux.*uxdx.*uxdy+4.*ux.^3.*uxdx.*uxdy+utdy.^2.*uy+(-1).*utdy.^2.*ux.^2.*uy+uxdy.^2.*uy+3.*ux.^2.*uxdy.^2.*uy+ux.*uxdydy.*uy+ux.^3.*uxdydy.*uy+(-1).*utdy.*uy.*uydt+utdy.*ux.^2.*uy.*uydt+uxdy.*uy.*uydx+3.*ux.^2.*uxdy.*uy.*uydx+ux.*uy.*uydxdy+ux.^3.*uy.*uydxdy+(-1).*ux.*uxdy.*uydy+ux.^3.*uxdy.*uydy+2.*ux.*uxdy.*uy.^2.*uydy+ux.*uydx.*uydy+ux.^3.*uydx.*uydy+(-2).*uy.*uydy.^2+2.*ux.^2.*uy.*uydy.^2+ut.*((-2).*utdt.*utdy.*((-1)+ux.^2)+(-1).*utdxdy.*(ux+ux.^3)+ux.*uxdtdy+ux.^3.*uxdtdy+(-1).*utdx.*uxdy+(-3).*utdx.*ux.^2.*uxdy+uxdt.*uxdy+3.*ux.^2.*uxdt.*uxdy+utdydy.*uy+(-1).*utdydy.*ux.^2.*uy+(-2).*utdy.*ux.*uxdy.*uy+2.*ux.*uxdy.*uy.*uydt+(-1).*uy.*uydtdy+ux.^2.*uy.*uydtdy+utdy.*uydy+(-1).*utdy.*ux.^2.*uydy+(-1).*uydt.*uydy+ux.^2.*uydt.*uydy)+(-1).*uydydy+(-1).*ux.^2.*uydydy+(-1).*uy.^2.*uydydy+ux.^2.*uy.^2.*uydydy);

    sUUd(:,:,2,3,1) = (1/2).*(uxdtdy+ux.^2.*uxdtdy+2.*ux.*uxdt.*uxdy+(-1).*utdtdt.*ux.*uy+ux.*uxdtdx.*uy+ux.^3.*uxdtdx.*uy+uxdt.*uxdx.*uy+3.*ux.^2.*uxdt.*uxdx.*uy+uxdtdy.*uy.^2+ux.^2.*uxdtdy.*uy.^2+2.*ux.*uxdt.*uxdy.*uy.^2+ux.*uxdx.*uydt+ux.^3.*uxdx.*uydt+2.*uxdy.*uy.*uydt+2.*ux.^2.*uxdy.*uy.*uydt+(-1).*ut.^2.*(utdtdt.*ux.*uy+utdt.*uxdt.*uy+utdt.*ux.*uydt)+(-1).*utdt.*(utdy.*ux.*(1+uy.^2)+uy.*(utdx.*(1+ux.^2)+(-1).*ux.*(ux.*uxdt+uy.*uydt)))+(-1).*ut.*(utdtdx.*uy+2.*utdt.^2.*ux.*uy+utdtdx.*ux.^2.*uy+2.*utdx.*ux.*uxdt.*uy+(-2).*ux.*uxdt.^2.*uy+(-1).*uxdtdt.*uy+(-1).*ux.^2.*uxdtdt.*uy+utdtdy.*ux.*(1+uy.^2)+utdx.*uydt+utdx.*ux.^2.*uydt+(-2).*uxdt.*uydt+(-1).*ux.^2.*uxdt.*uydt+(-1).*uxdt.*uy.^2.*uydt+(-2).*ux.*uy.*uydt.^2+utdy.*(uxdt+uxdt.*uy.^2+2.*ux.*uy.*uydt)+(-1).*ux.*uydtdt+(-1).*ux.*uy.^2.*uydtdt)+uydtdx+ux.^2.*uydtdx+uy.^2.*uydtdx+ux.^2.*uy.^2.*uydtdx+ux.*uy.*uydtdy+ux.*uy.^3.*uydtdy+2.*ux.*uxdt.*uydx+2.*ux.*uxdt.*uy.^2.*uydx+2.*uy.*uydt.*uydx+2.*ux.^2.*uy.*uydt.*uydx+uxdt.*uy.*uydy+uxdt.*uy.^3.*uydy+ux.*uydt.*uydy+3.*ux.*uy.^2.*uydt.*uydy);
    sUUd(:,:,2,3,2) = (1/2).*(uxdxdy+ux.^2.*uxdxdy+2.*ux.*uxdx.*uxdy+(-1).*utdtdx.*ux.*uy+(-1).*utdx.^2.*(1+ux.^2).*uy+(-1).*utdt.*uxdx.*uy+uxdx.^2.*uy+3.*ux.^2.*uxdx.^2.*uy+ux.*uxdxdx.*uy+ux.^3.*uxdxdx.*uy+uxdxdy.*uy.^2+ux.^2.*uxdxdy.*uy.^2+2.*ux.*uxdx.*uxdy.*uy.^2+utdx.*((1+ux.^2).*uxdt.*uy+(-1).*utdy.*ux.*(1+uy.^2)+ux.*(1+uy.^2).*uydt)+(-1).*utdt.*ux.*uydx+3.*ux.*uxdx.*uydx+ux.^3.*uxdx.*uydx+2.*uxdy.*uy.*uydx+2.*ux.^2.*uxdy.*uy.*uydx+2.*ux.*uxdx.*uy.^2.*uydx+2.*uy.*uydx.^2+2.*ux.^2.*uy.*uydx.^2+(-1).*ut.^2.*(utdtdx.*ux.*uy+utdt.*uxdx.*uy+utdt.*ux.*uydx)+ut.*((-1).*utdxdx.*uy+(-2).*utdt.*utdx.*ux.*uy+(-1).*utdxdx.*ux.^2.*uy+uxdtdx.*uy+ux.^2.*uxdtdx.*uy+(-2).*utdx.*ux.*uxdx.*uy+2.*ux.*uxdt.*uxdx.*uy+(-1).*utdxdy.*ux.*(1+uy.^2)+uxdx.*uydt+uxdx.*uy.^2.*uydt+ux.*uydtdx+ux.*uy.^2.*uydtdx+(-1).*utdx.*uydx+(-1).*utdx.*ux.^2.*uydx+uxdt.*uydx+ux.^2.*uxdt.*uydx+2.*ux.*uy.*uydt.*uydx+(-1).*utdy.*(uxdx+uxdx.*uy.^2+2.*ux.*uy.*uydx))+uydxdx+ux.^2.*uydxdx+uy.^2.*uydxdx+ux.^2.*uy.^2.*uydxdx+ux.*uy.*uydxdy+ux.*uy.^3.*uydxdy+uxdx.*uy.*uydy+uxdx.*uy.^3.*uydy+ux.*uydx.*uydy+3.*ux.*uy.^2.*uydx.*uydy);
    sUUd(:,:,2,3,3) = (1/2).*(2.*ux.*uxdy.^2+uxdydy+ux.^2.*uxdydy+(-1).*utdtdy.*ux.*uy+ux.*uxdxdy.*uy+ux.^3.*uxdxdy.*uy+(-1).*utdt.*uxdy.*uy+uxdx.*uxdy.*uy+3.*ux.^2.*uxdx.*uxdy.*uy+2.*ux.*uxdy.^2.*uy.^2+uxdydy.*uy.^2+ux.^2.*uxdydy.*uy.^2+(-1).*utdy.^2.*ux.*(1+uy.^2)+2.*ux.*uxdy.*uydx+2.*ux.*uxdy.*uy.^2.*uydx+uydxdy+ux.^2.*uydxdy+uy.^2.*uydxdy+ux.^2.*uy.^2.*uydxdy+(-1).*utdt.*ux.*uydy+ux.*uxdx.*uydy+ux.^3.*uxdx.*uydy+3.*uxdy.*uy.*uydy+2.*ux.^2.*uxdy.*uy.*uydy+uxdy.*uy.^3.*uydy+2.*uy.*uydx.*uydy+2.*ux.^2.*uy.*uydx.*uydy+ux.*uydy.^2+3.*ux.*uy.^2.*uydy.^2+(-1).*ut.^2.*(utdtdy.*ux.*uy+utdt.*uxdy.*uy+utdt.*ux.*uydy)+ut.*((-1).*utdxdy.*(1+ux.^2).*uy+uxdtdy.*uy+ux.^2.*uxdtdy.*uy+(-2).*utdx.*ux.*uxdy.*uy+2.*ux.*uxdt.*uxdy.*uy+(-1).*utdydy.*ux.*(1+uy.^2)+uxdy.*uydt+uxdy.*uy.^2.*uydt+ux.*uydtdy+ux.*uy.^2.*uydtdy+(-1).*utdx.*uydy+(-1).*utdx.*ux.^2.*uydy+uxdt.*uydy+ux.^2.*uxdt.*uydy+2.*ux.*uy.*uydt.*uydy)+(-1).*utdy.*(utdx.*(1+ux.^2).*uy+(-1).*uxdt.*uy+(-1).*ux.^2.*uxdt.*uy+(-1).*ux.*uydt+(-1).*ux.*uy.^2.*uydt+ut.*(uxdy.*(1+uy.^2)+2.*ux.*uy.*(utdt+uydy)))+ux.*uy.*uydydy+ux.*uy.^3.*uydydy);

    sUUd(:,:,3,3,1) = (1/2).*(utdt.*utdx.*ux+(-1).*utdt.*ux.*uxdt+(-1).*uxdtdx+(-1).*ux.^2.*uxdtdx+(-2).*ux.*uxdt.*uxdx+(-1).*utdt.*utdy.*uy+ux.*uxdtdy.*uy+uxdt.*uxdy.*uy+(-1).*utdt.*utdx.*ux.*uy.^2+utdt.*ux.*uxdt.*uy.^2+(-1).*uxdtdx.*uy.^2+ux.^2.*uxdtdx.*uy.^2+2.*ux.*uxdt.*uxdx.*uy.^2+(-1).*utdt.*utdy.*uy.^3+ux.*uxdtdy.*uy.^3+uxdt.*uxdy.*uy.^3+(-1).*utdtdt.*(1+uy.^2)+ux.*uxdy.*uydt+(-1).*utdt.*uy.*uydt+(-2).*uxdx.*uy.*uydt+2.*ux.^2.*uxdx.*uy.*uydt+3.*ux.*uxdy.*uy.^2.*uydt+utdt.*uy.^3.*uydt+ut.^2.*(utdtdt+(-1).*utdtdt.*uy.^2+(-2).*utdt.*uy.*uydt)+ut.*(utdtdx.*ux+utdx.*uxdt+(-1).*uxdt.^2+(-1).*ux.*uxdtdt+(-1).*utdtdy.*uy+(-1).*utdtdx.*ux.*uy.^2+(-1).*utdx.*uxdt.*uy.^2+uxdt.^2.*uy.^2+ux.*uxdtdt.*uy.^2+(-1).*utdtdy.*uy.^3+(-2).*utdt.^2.*((-1)+uy.^2)+(-1).*utdy.*uydt+(-2).*utdx.*ux.*uy.*uydt+2.*ux.*uxdt.*uy.*uydt+(-3).*utdy.*uy.^2.*uydt+uydt.^2+3.*uy.^2.*uydt.^2+uy.*uydtdt+uy.^3.*uydtdt)+ux.*uy.*uydtdx+ux.*uy.^3.*uydtdx+uydtdy+2.*uy.^2.*uydtdy+uy.^4.*uydtdy+uxdt.*uy.*uydx+uxdt.*uy.^3.*uydx+ux.*uydt.*uydx+3.*ux.*uy.^2.*uydt.*uydx+4.*uy.*uydt.*uydy+4.*uy.^3.*uydt.*uydy);
    sUUd(:,:,3,3,2) = (1/2).*(utdx.^2.*ux+(-1).*utdx.*ux.*uxdt+(-2).*ux.*uxdx.^2+(-1).*uxdxdx+(-1).*ux.^2.*uxdxdx+(-1).*utdx.*utdy.*uy+ux.*uxdxdy.*uy+uxdx.*uxdy.*uy+(-1).*utdx.^2.*ux.*uy.^2+utdx.*ux.*uxdt.*uy.^2+2.*ux.*uxdx.^2.*uy.^2+(-1).*uxdxdx.*uy.^2+ux.^2.*uxdxdx.*uy.^2+(-1).*utdx.*utdy.*uy.^3+ux.*uxdxdy.*uy.^3+uxdx.*uxdy.*uy.^3+(-1).*utdtdx.*(1+uy.^2+ut.^2.*((-1)+uy.^2))+utdx.*uy.*uydt+utdx.*uy.^3.*uydt+ux.*uxdy.*uydx+(-2).*utdt.*uy.*uydx+(-2).*ut.^2.*utdt.*uy.*uydx+(-1).*uxdx.*uy.*uydx+2.*ux.^2.*uxdx.*uy.*uydx+3.*ux.*uxdy.*uy.^2.*uydx+uxdx.*uy.^3.*uydx+ux.*uydx.^2+3.*ux.*uy.^2.*uydx.^2+ut.*(utdxdx.*ux+(-1).*ux.*uxdtdx+utdx.*uxdx+(-1).*uxdt.*uxdx+(-1).*utdxdy.*uy+(-1).*utdxdx.*ux.*uy.^2+ux.*uxdtdx.*uy.^2+(-1).*utdx.*uxdx.*uy.^2+uxdt.*uxdx.*uy.^2+(-1).*utdxdy.*uy.^3+(-2).*utdt.*utdx.*((-1)+uy.^2)+uy.*uydtdx+uy.^3.*uydtdx+(-1).*utdy.*uydx+(-2).*utdx.*ux.*uy.*uydx+2.*ux.*uxdt.*uy.*uydx+(-3).*utdy.*uy.^2.*uydx+uydt.*uydx+3.*uy.^2.*uydt.*uydx)+ux.*uy.*uydxdx+ux.*uy.^3.*uydxdx+uydxdy+2.*uy.^2.*uydxdy+uy.^4.*uydxdy+4.*uy.*uydx.*uydy+4.*uy.^3.*uydx.*uydy);
    sUUd(:,:,3,3,3) = (1/2).*(utdx.*utdy.*ux+(-1).*utdy.*ux.*uxdt+(-1).*uxdxdy+(-1).*ux.^2.*uxdxdy+(-2).*ux.*uxdx.*uxdy+(-1).*utdy.^2.*uy+uxdy.^2.*uy+ux.*uxdydy.*uy+(-1).*utdx.*utdy.*ux.*uy.^2+utdy.*ux.*uxdt.*uy.^2+(-1).*uxdxdy.*uy.^2+ux.^2.*uxdxdy.*uy.^2+2.*ux.*uxdx.*uxdy.*uy.^2+(-1).*utdy.^2.*uy.^3+uxdy.^2.*uy.^3+ux.*uxdydy.*uy.^3+(-1).*utdtdy.*(1+uy.^2+ut.^2.*((-1)+uy.^2))+utdy.*uy.*uydt+utdy.*uy.^3.*uydt+uxdy.*uy.*uydx+uxdy.*uy.^3.*uydx+ux.*uy.*uydxdy+ux.*uy.^3.*uydxdy+ux.*uxdy.*uydy+(-2).*utdt.*uy.*uydy+(-2).*ut.^2.*utdt.*uy.*uydy+(-2).*uxdx.*uy.*uydy+2.*ux.^2.*uxdx.*uy.*uydy+3.*ux.*uxdy.*uy.^2.*uydy+ux.*uydx.*uydy+3.*ux.*uy.^2.*uydx.*uydy+4.*uy.*uydy.^2+4.*uy.^3.*uydy.^2+ut.*(utdxdy.*ux+(-1).*ux.*uxdtdy+utdx.*uxdy+(-1).*uxdt.*uxdy+(-1).*utdydy.*uy+(-1).*utdxdy.*ux.*uy.^2+ux.*uxdtdy.*uy.^2+(-1).*utdx.*uxdy.*uy.^2+uxdt.*uxdy.*uy.^2+(-1).*utdydy.*uy.^3+(-2).*utdt.*utdy.*((-1)+uy.^2)+uy.*uydtdy+uy.^3.*uydtdy+(-1).*utdy.*uydy+(-2).*utdx.*ux.*uy.*uydy+2.*ux.*uxdt.*uy.*uydy+(-3).*utdy.*uy.^2.*uydy+uydt.*uydy+3.*uy.^2.*uydt.*uydy)+uydydy+2.*uy.^2.*uydydy+uy.^4.*uydydy);


    % wUU, antisymmetric
    wUU = zeros(params.Nx, params.Ny, 3, 3);
    xyzero = zeros(params.Nx, params.Ny);
    wUU(:,:,1,2) = (1/2).*(utdx.*((-1)+ut.^2+(-1).*ux.^2)+((-1)+ut.^2+(-1).*ux.^2).*uxdt+(-1).*uy.*(utdy.*ux+(-1).*ut.*uxdy+ux.*uydt+ut.*uydx));
    wUU(:,:,1,3) = (1/2).*((-1).*utdx.*ux.*uy+(-1).*ux.*uxdt.*uy+utdy.*((-1)+ut.^2+(-1).*uy.^2)+(-1).*uydt+ut.^2.*uydt+(-1).*uy.^2.*uydt+ut.*ux.*((-1).*uxdy+uydx));
    wUU(:,:,2,3) = (1/2).*(ut.*(utdy.*ux+(-1).*utdx.*uy+(-1).*uxdt.*uy+ux.*uydt)+(-1).*(1+ux.^2+uy.^2).*(uxdy+(-1).*uydx));


    PiUUdd = zeros(params.Nx, params.Nx, 3, 3, 3, 3); % fml
    PiUUdd(:,:,1,1,1,:) = reshape([(1/2).*((-1)+ut.^2).^2,(1/2).*(ut.*ux+(-1).*ut.^3.*ux),(1/2).*(ut.*uy+(-1).*ut.^3.*uy)], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,1,1,2,:) = reshape([(1/2).*(ut.*ux+(-1).*ut.^3.*ux),(1/2).*(1+ux.^2+ut.^2.*((-1)+ux.^2)),(1/2).*(1+ut.^2).*ux.*uy], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,1,1,3,:) = reshape([(1/2).*(ut.*uy+(-1).*ut.^3.*uy),(1/2).*(1+ut.^2).*ux.*uy,(1/2).*(1+uy.^2+ut.^2.*((-1)+uy.^2))], Nx, Ny, 1, 1, 1, 3);

    PiUUdd(:,:,1,2,1,:) = reshape([(1/2).*((-1)+ut).*ut.*(1+ut).*ux,(-1/2).*((-1)+ut).*(1+ut).*(1+ux.^2),(-1/2).*((-1)+ut).*(1+ut).*ux.*uy], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,1,2,2,:) = reshape([(-1/2).*((-1)+ut).*(1+ut).*(1+ux.^2),(1/2).*ut.*ux.*(1+ux.^2),(1/2).*ut.*(1+ux.^2).*uy], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,1,2,3,:) = reshape([(-1/2).*((-1)+ut).*(1+ut).*ux.*uy,(1/2).*ut.*(1+ux.^2).*uy,(1/2).*ut.*ux.*((-1)+uy.^2)], Nx, Ny, 1, 1, 1, 3);

    PiUUdd(:,:,1,3,1,:) = reshape([(1/2).*((-1)+ut).*ut.*(1+ut).*uy,(-1/2).*((-1)+ut).*(1+ut).*ux.*uy,(-1/2).*((-1)+ut).*(1+ut).*(1+uy.^2)], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,1,3,2,:) = reshape([(-1/2).*((-1)+ut).*(1+ut).*ux.*uy,(1/2).*ut.*((-1)+ux.^2).*uy,(1/2).*ut.*ux.*(1+uy.^2)], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,1,3,3,:) = reshape([(-1/2).*((-1)+ut).*(1+ut).*(1+uy.^2),(1/2).*ut.*ux.*(1+uy.^2),(1/2).*ut.*uy.*(1+uy.^2)], Nx, Ny, 1, 1, 1, 3);

    PiUUdd(:,:,2,2,1,:) = reshape([(1/2).*(1+ux.^2+ut.^2.*((-1)+ux.^2)),(-1/2).*ut.*ux.*(1+ux.^2),(-1/2).*ut.*((-1)+ux.^2).*uy], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,2,2,2,:) = reshape([(-1/2).*ut.*ux.*(1+ux.^2),(1/2).*(1+ux.^2).^2,(1/2).*ux.*(1+ux.^2).*uy], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,2,2,3,:) = reshape([(-1/2).*ut.*((-1)+ux.^2).*uy,(1/2).*ux.*(1+ux.^2).*uy,(1/2).*((-1)+(-1).*uy.^2+ux.^2.*((-1)+uy.^2))], Nx, Ny, 1, 1, 1, 3);

    PiUUdd(:,:,2,3,1,:) = reshape([(1/2).*(1+ut.^2).*ux.*uy,(-1/2).*ut.*(1+ux.^2).*uy,(-1/2).*ut.*ux.*(1+uy.^2)], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,2,3,2,:) = reshape([(-1/2).*ut.*(1+ux.^2).*uy,(1/2).*ux.*(1+ux.^2).*uy,(1/2).*(1+ux.^2).*(1+uy.^2)], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,2,3,3,:) = reshape([(-1/2).*ut.*ux.*(1+uy.^2),(1/2).*(1+ux.^2).*(1+uy.^2),(1/2).*ux.*uy.*(1+uy.^2)], Nx, Ny, 1, 1, 1, 3);

    PiUUdd(:,:,3,3,1,:) = reshape([(1/2).*(1+uy.^2+ut.^2.*((-1)+uy.^2)),(-1/2).*ut.*ux.*((-1)+uy.^2),(-1/2).*ut.*uy.*(1+uy.^2)], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,3,3,2,:) = reshape([(-1/2).*ut.*ux.*((-1)+uy.^2),(1/2).*((-1)+(-1).*uy.^2+ux.^2.*((-1)+uy.^2)),(1/2).*ux.*uy.*(1+uy.^2)], Nx, Ny, 1, 1, 1, 3);
    PiUUdd(:,:,3,3,3,:) = reshape([(-1/2).*ut.*uy.*(1+uy.^2),(1/2).*ux.*uy.*(1+uy.^2),(1/2).*(1+uy.^2).^2], Nx, Ny, 1, 1, 1, 3);

    % fill the symmetric part
    for i=1:3
      for j=i+1:3
        sUU(:,:,j,i)    =  sUU(:,:,i,j);
        sUUd(:,:,j,i,:) = sUUd(:,:,i,j,:);
        PiUUdd(:,:,j,i,:,:) = PiUUdd(:,:,i,j,:,:);
        wUU(:,:,j,i) = -wUU(:,:,i,j);
      end
    end

    % d_a sig^{mu nu} u^a
    dersuUU = zeros(params.Nx, params.Nx, 3, 3);
    for i=1:3 
      for xi=1:params.Nx
        for yi=1:params.Ny
          %size(sUUd)
          %size(uU)
          %size(sUUd(xi,yi,:,:,i) .* uU(xi,yi,i))
          dersuUU(xi,yi,:,:) = dersuUU(xi,yi,:,:) + reshape(sUUd(xi,yi,:,:,i) .* uU(xi,yi,i), 1, 1, 3, 3);
        end
      end
    end

    divu = utdt+uxdx+uydy;

    S1UU = zeros(params.Nx, params.Ny, 3, 3);
    for xi=1:Nx
      for yi=1:Ny
        S1UU(xi,yi,:,:) = dersuUU(xi,yi,:,:) + 1/2 .* sUU(xi,yi,:,:) .* divu(xi,yi);
      end
    end
    for xi=1:params.Nx
      for yi=1:params.Ny
        tmp = zeros(1,1,3,3);
        for i=1:3
          for j=1:3
            tmp = tmp + PiUUdd(xi,yi,:,:,i,j) .* S1UU(xi,yi,i,j);
          end
        end
        S1UU(xi,yi,:,:) = tmp;
      end
    end

    sdU = sUU;
    sdU(:,:,1,:) = -sUU(:,:,1,:);

    S2UU = zeros(Nx, Ny, 3, 3);
    for xi=1:Nx
      for yi=1:Ny
        for i=1:3
          for c=1:3
            S2UU(xi,yi,i,:) = S2UU(xi,yi,i,:) + reshape(sdU(xi,yi,c,i) .* wUU(xi,yi,:,c), 1, 1, 1, 3);
          end 
        end
      end
    end

    for xi=1:params.Nx
      for yi=1:params.Ny
        tmp = zeros(1,1,3,3);
        for i=1:3
          for j=1:3
            tmp = tmp + PiUUdd(xi,yi,:,:,i,j) .* S2UU(xi,yi,i,j);
          end
        end
        S2UU(xi,yi,:,:) = tmp;
      end
    end

    T0UU = zeros(Nx,Ny,3,3);
    T1UU = zeros(Nx,Ny,3,3);
    T2UU = zeros(Nx,Ny,3,3);

    T0UU(:,:,1,:) = reshape([(1/2).*T.^3.*((-1)+3.*ut.^2),(3/2).*T.^3.*ut.*ux,(3/2).*T.^3.*ut.*uy], Nx, Ny, 1, 3);
    T0UU(:,:,2,:) = reshape([(3/2).*T.^3.*ut.*ux,(1/2).*T.^3.*(1+3.*ux.^2),(3/2).*T.^3.*ux.*uy], Nx, Ny, 1, 3);
    T0UU(:,:,3,:) = reshape([(3/2).*T.^3.*ut.*uy,(3/2).*T.^3.*ux.*uy,(1/2).*T.^3.*(1+3.*uy.^2)], Nx, Ny, 1, 3);

    for xi=1:Nx
      for yi=1:Ny
        T1UU(xi,yi,:,:) = T1UU(xi,yi,:,:) - T(xi,yi).^2 .* sUU(xi,yi,:,:);
      end
    end

    for xi=1:params.Nx
      for yi=1:params.Ny
        T2UU(xi,yi,:,:) = T2UU(xi,yi,:,:) + T(xi,yi).*( (params.A+1).*S1UU(xi,yi,:,:) + 2.*params.A.*S2UU(xi,yi,:,:) );
      end
    end




end

